Damage spreading in small world Ising models 
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We study damage-spreading in the ferromagnetic Ising model on small world networks 
using Monte Carlo simulation with Glauber dynamics. The damage spreading temperature 
T d is determined as a function of rewiring probability p for small world networks obtained 
by rewiring the 2D square and 3D cubic lattices. We find that the damage for different 
values of p collapse onto master curves when plotted against a rescaled temperature and 
that the distance between T d and the critical temperature T c increases with p. We argue 
that when using the Ising model to study social systems, it is necessary to place the spins 
on a small world network rather than on a regular lattice. 

PACS numbers: 75.10.Hk (Classical spin models), 75.40.Mg (Numerical simulation 
studies), 75.10.Nr (Spin-glass and other random models) 



I. INTRODUCTION 

The Ising model is one of the most Important 
models of statistical mechanics. It and its gen- 
eralisations have been used to model a variety of 
natural phenomena, ranging from biology to com- 
puter science and social science (e.g., [[lj, ^ |3|. |5|j). 
For instance, many social systems can be modelled 
by letting spin up/down denote different opinions 
or preferences. In such models, a ferromagnetic 
interaction is interpreted as two people who prefer 
to agree, while an antiferromagnetic interactions 
means that they want to disagree. A magnetic field 
adds a bias than can be interpreted as "prejudices" 
or "stubbornness", while the randomness induced 
by a finite temperature can be seen as a "free will". 

Damage spreading is a tool for studying the in- 
fluence of perturbations on the equilibrium state 
of a system. It has been used to determine some 
properties of the energy landscape for disordered 
spin systems [pi , and also has great uses for play- 
ing "what if '-type scenarios in models of complex 
systems. For a voter model, for instance, damage 
spreading studies how much influence a (small) 
set of voters can have over the final outcome of 
the election. Damage spreading was first used by 
Kauffman |J6|] as a tool for studying biologically mo- 
tivated dynamical systems, but has since found 



widespread use also in physics (e.g., (?]]). 

Damage spreading works by duplicating an equi- 
librium spin configuration of a system and chang- 
ing a fraction d of the spins. Both systems are 
then subjected to the same thermal noise and the 
distance between them is calculated. In Monte 
Carlo simulations, both systems are simulated si- 
multaneously: the same spin is selected for spin- 
flip in both systems, and the same random number 
("thermal noise") is used to determine whether an 
energy-raising flip should be performed. 

After equilibrating both systems, the Hamming 
distance (the number of different spins) between 
the spin configurations S a and S 13 
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(where 8 is the Kronecker delta function) is mea- 
sured. The Hamming distance can also be ex- 
pressed in terms of the Parisi overlap [pi 
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Most of the work on both spin models and dam- 
age spreading place the spins either on a finite- 
dimensional lattice or on a random grap h. Here 
we instead use small world graphs p|, [XOQ to study 
the ferromagnetic Ising model on graphs interpo- 
lating between 2 and 3-dimensional simple cubic 
lattices and random graphs with the same connec- 



tivity. The Hamiltonian of our model is 



H — 2_^ JijSiSj 

i<j 



(3) 



where J tj is 1 if and only if there is an edge between 
spins i and j and otherwise. 

Small world graphs are intermediates between 
a regular lattice and a random graph; they have 
previously been used to study, e.g., computation, 
diffusion, and spreading of diseases. The origi- 
nal motivation for studying small worlds is that 
they possess both small diameters (like a random 
graph [11|) and a high degree of clustering (like 
a regular lattice). For examples of real world net- 
works with small world characteristics and reviews 
of previous work, see, e.g., JTo|, |l2|, 13, 14|. 



The small world is constructed by considering 
in turn all the edges (i,j) of a lattice and with 
some probability p replacing it with a random edge 
(i,k). The rewiring parameter p thus determines 
how many of the links are removed and can be 
used to interpolate between the regular lattice and 
a random graph. Note that the small world for 
p = 1 differs slightly from a random graph, since 
all nodes are guaranteed to have a local connectiv- 
ity of at least z/2 where z is the connectivity of the 
regular lattice. The distribution of connectivities 
is more broad for the small world with p = 1. We 
chose to use the small world model where links are 
rewired and not the one where they are added be- 
cause we wanted to keep the average connectivity 
of the graphs the same for all p. 

The use of small world graphs to study physi- 
cal models has so far been limited. Barrat and 
Weigt 1 15 1 and Gitterman [ |l6| ] have used them 
to studythe crossover from ID to mean-field be- 
haviour for the ferromagnetic Ising model, finding 
a disorder -order transition at a finite temperature 
T c (p) for any p > 0, provided that the system size is 
large enough. 

Most of the work on small world networks has 
started by rewiring a one-dimensional ring lattice, 
but here we instead use the 2D square and 3D 
simple cubic lattices. One reason for doing this 
is that while the ID Ising model is trivial and dis- 
ordered for all finite temperatures, the 2D and 3D 
versions are ordered below a critical temperature 
T c . The 2D model can be solved exactly, while for 
4D and higher-dimensions, mean field theory ex- 
plains the phase transition (see, e.g., [17]). An 



important concept in the study of phase transi- 
tions and critical phenomena is that of universality 
class. Models displaying the same behaviour close 
to T c are said to be in the same universality class, 
and it turns out that there are many fewer uni- 
versality classes than models. Putting spin models 
on small world graphs provides an opportunity to 



study the crossover from a finite-dimensional uni- 
versality class to mean field behaviour. Here we 
restrict ourselves to determining T c , but it would 
also be interesting to see how the critical expo- 
nents change as p is increased. 

It should be noted that the small world networks 
used here differ from those obtained by rewiring 
a ring lattice in one respect: their clustering co- 
efficient does not display the same threshold be- 
haviour as a function of p: it starts at for p = 
(since the regular lattices used are bipartite) and 
then grows to the random graph value. The graphs 
used here are however still clustered in the sense 
that if j and k are neighbours of i, then there 
is a short path between them that does not pass 
through i. 

While the emphasis in the present work is on 
the damage spreading behaviour of the model, we 
also determined the critical temperature T c for the 
order -disorder transition. This was done primarily 
in order to compare it with the damage spreading 
temperature Td\ the numerical accuracy of T c is 
smaller than that for T d . 

The Monte Carlo method used was the standard 



single spin-flip Metropolis | IS] algorithm. In each 
time-step, N spin flips are attempted. For each 
flip-attempt, a spin is randomly selected and the 
energy-change AH if it is flipped is calculated. If 
the change in energy is negative, the spin is al- 
ways flipped, otherwise it is flipped with probabil- 
ity e - AH / T where T is the temperature. We also did 
some runs using different MC procedures (heat- 
bath algorithm, spin-exchange, using an ordered 
update instead of a random) . We found that using 
the heat-bath algorithm caused the damage to heal 
at temperatures close to and above T c , while for the 
spin-exchange dynamics with the Metropolis algo- 
rithm the damage spreads for all temperatures. 
Updating the spins in order instead of randomly 
gives a smaller damage for all temperatures. These 
results agree with the results of Vojta [IS, 20, 21 1 
for the standard Ising model. 

In most of the simulations, we used the Mitchell- 
Moore additive random number generator (see, 
e.g., |22| for a description). We also did some runs 
with the standard C library's drand4 8 () genera- 
tor and found the same behaviour. All simula- 
tions were averaged over JV| different rewiring pro- 
cedures, and for each small world graph an av- 
erage over N r independent Monte Carlo runs was 
performed. Typical values were JVj = N r = 10, but 
this was varied for some runs in order to check 
self-averaging. No significant differences in be- 
haviour was found. 

Our simulation procedure was simple. After 
equilibrating the system (using simulated anneal- 
ing), a copy is made and d N spins in it are flipped. 
Both systems are then simulated using the same 
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FIG. 1: This figure shows the damage as a function of 
temperature for small world graphs obtained by rewiring 
a 100 x 100 2D square lattice with (from left to right) p — 0, 
0.01, 0.05, 0.1, 0.2, 0.4, 0.6, and 1. For each p, an av- 
erage over 10 graphs and 10 restarts per graph was per- 
formed. The location of Td shifts to higher temperatures 
as p is increased, and the slope of d(T) decreases. 



random numbers to determine which spin to se- 
lect and whether or not to flip it. After equilib- 
rium has been reached again, we start measuring 
the damage as well as other quantities such as the 
magnetisation and energy and their standard devi- 
ations. We used d = 0.01 in all of the simulations 
presented here; none of the results presented are 
sensitive to the exact value of d a . In order to check 
the dependence on initial conditions, we also per- 
formed some runs damaging a non-equilibrated 
system; these gave the same results. 

Figure [I] shows the end-damage as a function of 
temperature for p ranging from to 1 . The rewired 
lattice in this figure is the 2D square with TV = 10 4 
spins. We tested some different system sizes and 
found that this seems to be a large enough num- 
ber of spins that finite-size effects are minimised. 
The data was averaged over Ni = 10 graphs and 
for each graph the Monte Carlo simulation was 
restarted N r = 10 times in order to improve nu- 
merical accuracy. Error bars for the damage in 
this and the following figures were determined to 
be at most on the order of 0.01 and in almost all 
cases considerably smaller. Note though that the 
errors increase with p, as should be expected since 
the averaging becomes more important for large p. 

Figure g] shows the corresponding data for the 
3D lattice. The system size here is N = 8000 and 
N t = N r = 10 as for the 2D data. 

We can define a damage spreading temperature 
T d {e) as the lowest temperature for which the dam- 
age d is larger than some (small) e, 



T d (e) = min{T : d(T) > e}. 



(4) 



FIG. 2: Here we show the damage as a function of tem- 
perature for small world graphs obtained by rewiring a 
20 x 20 x 20 3D cubic lattice with (from left to right) p = 0, 
0.01, 0.05, 0.1, 0.2, 0.4, 0.6, and 1. For each p, an 
average over 10 graphs and 10 restarts per graph was 
performed. As in the 2d case, the location of T d shifts to 
higher temperatures as p is increased, and the slope of 
d(T) decreases. 
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In the limit as e — > 0, our T d (e) converges to the 
standard T d which is defined as the lowest temper- 
ature for which the damage is non-zero. We use a 



FIG. 3: For the same data as in figure m this figure 
shows the p-dependence of T c (squares) and T d for some 
different e. Note the logarithmic scale of the p-axis in 
this plot. It is clear that Td is independent of e for small 
enough e's. 



non-zero e in equation ^ when determining T d from 
our data becuase using a e smaller than the error - 
bar for the damage would lead to noise in T d . Fig- 
ures [I] and ^ show clearly that T d increases with 
p, as is to be expected. In order to quantify this, 
figure m compares T d to the order -disorder tran- 
sition temperature T c for the 2D data. The figure 
shows T d for e = 10~ 4 , 10~ 3 , 10~ 2 , and 10 _1 ; it is 
clear that the definition of T d is independent of e 
for small enough e's. The temperature where the 
damage attained its maximum value of 0.5 seems 
to approach T c ; this is in agreement with previous 
work [23]. The critical temperature T c was deter- 
mined as the temperature at which the Binder's 
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FIG. 4: Here we show T c (squares) and T d as a function 
of p for some different e for the 3d case. Here, too, the 
values for T d are independent of the exact value of e, 
provided that it is small enough. 



cumulant 



(to 2 ) 2 



(5) 



curves for large system sizes cross. For the 2D 
lattice, the largest system simulated consisted of 
10 4 spins, while in the 3D case shown in figure § 
below, system sizes up to 21 3 = 9261 were used 
to determine T c . The error bars for T c are larger 
than for T d \ note that the mean-field value for (reg- 
ular) random graphs with coordination number z 
is T c — z. The value of T d obtained for p = 1 here 
is in reasonable agreement to the one for normal 
random graphs. 



2D 
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0.01 


0.05 


0.1 


0.2 


0.4 
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T d 


2.24 


2.28 


2.34 


2.40 


2.49 


2.60 


2.70 


2.83 


3D 


P 





0.01 


0.05 


0.1 
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1 




T d 


4.08 


4.08 


4.10 


4.14 


4.19 


4.29 


4.36 


4.48 



TABLE I: T d for the small world starting from a 2D and 
3D lattices. 



Table § shows the values for T d for different p for 
small worlds obtained by rewiring the 2D square 
and 3D cubic lattices. For p = 0, we get values in 
agreement with those reported in the literature! 20, 
24, f25h. 



Scaling plots are used to combine data from runs 
with different values of some parameter into one 
curve. In our case, we can make the data for dif- 
ferent p fall onto the same curve by plotting the 
damage as a function of a rescaled temperature 



T 



T-T d 

A(p) 



(6) 



FIG. 5: Same data as in figure H plotted as a function 
of a rescaled temperature. By plotting the damage as a 
function of a reduced temperature T — (T — T d )/A(p), it 
is possible to get collapse for all p except p — (shown as 
small squares in the figure), which does not follow the 
same functional form as the other curves. 



Our scaling ansatz is that the damage can be writ- 
ten as 



D(T,p) = /( 



T-T d {p) , 
A(p) ' 



(7) 



for some / which is independent of p. In equa- 
tion |], A(p) is determined by the inverse of the 
rate at which the damage develops for different p 



1 df 



dD 

lT {T - Td) A(p) d f 



~(T = 0). 



(8) 



A is an increasing function of p; physically it tells 
us how much more we must increase the temper- 
ature in order to get the same increase in damage 
for different p: 



AT oc A(p)AD. 



(9) 



The values for A(p) determined from the data in 
figures [l] and are shown in table O. We found a 
reasonable scaling A(p) ~ p a with a rs 0.35 for the 
2D data and a w 0.2 for the 3D data. The function 
/ turns out to be linear. 



2D 
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0.01 


0.05 


0.1 
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A(p) 


0.07 


0.1 


0.18 


0.25 


0.3 


0.38 


0.43 


0.5 


3D 


P 





0.01 


0.05 


0.1 


0.2 


0.4 
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1 




A(p) 


0.28 


0.3 


0.36 


0.41 


0.5 


0.55 


0.61 


0.7 



TABLE II: A(p) for the 2D and 3D rewired lattices. 

Figure ||| plots the damage as a function of f for 
the 2D case. A very good collapse is obtained for 
all p > 0. The data for p = can not be made to 
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FIG. 6: In contrast to the 2D case, by plotting the dam- 
age as a function of a reduced temperature (T—Td)/A(p), 
it is possible to get collapse for all p for the 3D data. 



fall onto the same curve. Note that the distance 
between the master curve and the p = data is 
larger than the estimated error bars. 

Figure ^ shows that, in contrast to the 2D case, 
the 3D data do collapse onto one curve for all p, 
including the p = (i.e., simple cubic lattice) case. 

This shows some qualitative differences between 
the 2 and 3-dimensional lattices. The way that 
damage spreads in the model can be seen as a form 
of generalised random walk; we speculate that the 
difference between the 2D p = data and the other 
data might be related to the differences (in, e.g., 
return time) between random walks on 2D and 



3D /random lattices [26|. 



We also studied the approach to equilibrium of 
the damaged system. Figure M below shows the re- 
laxation of the damage as a function of the number 
of complete Monte Carlo sweeps after the damage 
is introduced. The figure shows data for 2D model 
with p — 0.4; the relaxation behaviour for other val- 
ues of p as well as for the 3D case is similar. It is 
clearly seen that there is a power -law for a short 
interval above Td. 

The data can be very approximately fitted to a 
form d(t) ~ t a with a « 1.5 ± 0.1 for T considerably 
larger than T d and for all p > 0. The exponent for 
p = is significantly different, a w 1.1. 

In conclusion, we found that the damage for 
different small worlds fall onto a universal curve 
when plotted as a function of a rescaled tempera- 



ture. The distance between T d and T c increases as 
a function of rewiring probability p, i.e., the range 
in temperature where the model is ordered but 
small perturbations are important increases. This 
is important for models of social systems, where 
we can interpret the temperature as a form of (ran- 
dom) "free will". 

We believe that putting spin models on small 
world graphs provides an ideal method not only of 




(MC sweeps per spin) 



FIG. 7: This figure shows the time-dependence of the 
damage for the 2D model with p = 0.4 and T = 2.2, 2.4, 
2.6, 2.8, and 3.0. The relaxation is exponential below Td, 
and displays a power -law for a short interval for T > T d . 
The damage spreading transition takes place at T d fa 2.6. 



studying social models more realistically but also 
of testing hypotheses regarding spin models. For 
instance, it is an interesting open question how 
to accurately describe the ground state and low- 
lying excitations of the 3D ±J spin glass model. 
By putting this model on a small world graph and 
studying the crossover to the p = 1 mean-field be- 
haviour, it might be possible to learn more about 
this. 
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